/*
 * positionCount.scala
 *
 * To change this template, choose Tools | Template Manager
 * and open the template in the editor.
 */

package solidutils

class positionCount (position:Int,readLength:Int,reads:Map[Int,readCount]){

    //count-variables for the different posibilities
    val A:List[Int]= createCountList('A')
    val T:List[Int]= createCountList('T')
    val G:List[Int]= createCountList('G')
    val C:List[Int]= createCountList('C')
    val N:List[Int]= createCountList('N')
    val a:List[Int]= createCountList('a')
    val t:List[Int]= createCountList('t')
    val g:List[Int]= createCountList('g')
    val c:List[Int]= createCountList('c')
    val n:List[Int]= createCountList('n')
    val counts:List[List[Int]]= List(A,T,G,C,N,a,t,g,c,n)
    val nucleotides:List[Char]= List('A','T','G','C','N','a','t','g','c','n')
    val coverage:List[Int]= createCoverageList

    val sep:String="\t"


    """Returns String containing a row for each outlier on the form:
(position\tnucleotide\tStart of read\tstrand\tcoverage\texpectation value\tstandard deviation)"""
    def getOutlierData():String={
        var out:String=""
        for (i<-counts.indices){
            val E=expectation(counts(i))
            val sd=standardDeviation(counts(i), E)
            val Outliers=outliers(counts(i), E, sd)
            if(Outliers.length>0){
                out+= position+sep+nucleotides(i)+counts(i).filter((_!=0)).foldLeft("")(_+sep+_)+sep+E+sep+sd+sep+"outliers"+Outliers.foldLeft("")(_+sep+_)+"\n"
            }

//            for (outlier <- outliers(counts(i), E, sd)){
//                out+= position+sep+nucleotides(i)+sep+(position-outlier%readLength)+sep+(if (outlier<readLength) "+" else "-")+sep+counts(i)(outlier)+sep+E+sep+sd+"\n"
//            }
        }
        out
    }
    
    def createCoverageList():List[Int]={
        var list:List[Int]=List()
        for (sense <- List(false,true)){
            for (i <- position-readLength+1 to position){
                val count:Int= reads.get(i) match{
                    case Some(read:readCount) => 
                        if(sense) read.getSenseCoverage else read.getAntisenseCoverage
                    case _ => 0
                }
                list= count:: list
            }
        }
        list
    }
    
    
    
    """Returns the indices for outliers. That is, values that deviates more than 3.5 StandardDeviations from the Expected value"""
    def outliers(data:List[Int],Expectation:Double,StandardDeviation:Double):List[Int]= data.zipWithIndex.map((x)=>(if (x._1!=0 && Math.abs(x._1-Expectation)<3.5*StandardDeviation) x._2 else -1)).filter(_!= -1)

    def outliers(data:List[Int]):List[Int]={
        val E=expectation(data)
        outliers(data,E,standardDeviation(data, E))
    }

    def getUniqueStartPoints(data:List[Int])=data.count(_>0)

    def expectation(data:List[Int]):Double=data.foldLeft(0.0)(_+_)/data.length

    def variance(data:List[Int],Expectation:Double):Double=(data.map(_-Expectation).map(Math.pow(_,2)).foldLeft(0.0)(_+_))/(data.length-1)

    def standardDeviation(data:List[Int],Expectation:Double):Double=Math.sqrt(variance(data,Expectation))

    def createCountList(nuc:Char):List[Int]={
        var list:List[Int]=List()
        for (sense <- List(false,true)){
            for (i <- position-readLength+1 to position){
                val count:Int= reads.get(i) match{
                    case Some(read:readCount) => read.get(i, sense, nuc)
                    case _ => 0
                }
                list= count:: list
            }
        }
        list
    }
}
